

beta=betavec(Qbeta)

global gammagrid elementsize lowerbound n agent P beta S state agent Y scalefactor rho ewdev statedev12 statedev13 statedev23
agent=2;
pen=1;
elementsize=0.1;
lowerbound=0.2;
gammagrid=[];
gammagrid(:,1)=[lowerbound:elementsize:1-lowerbound];
gammagrid(:,2)=1-gammagrid(:,1);

gammagrid(:,3)=1;
gammagrid(:,4)=gammagrid(:,2)./gammagrid(:,1);
n=length(gammagrid);



penalty=0+0.01*(pen-1);

income=incomevals;
S=gridmake(income,income);
SS=gridmake([1:length(income(:,1)) ]', [1:length(income(:,1)) ]');

state=length(S);
Y=sum(S,2);	
% possible aggregate income outcome
P=zeros(state,state);
for k=1:state
    for j=1:state
        P(j,k)=Q1(SS(j,1),SS(k,1))*Q1(SS(j,2),SS(k,2));
    end
end

agent=2;  

% Compute set of autarky values
caut=zeros(state,agent);
for k=1:agent
caut(:,k)=S(:,k);
end
waut=zeros(state,agent);
wautnew=zeros(state,agent),
dwaut=1;
crit=0.0000000001;
while dwaut>crit
    
    wautnew(:,1)=utility(rho,caut(:,1))+beta*P*waut(:,1);
    wautnew(:,2)=utility(rho,caut(:,2))+beta*P*waut(:,2);
    dwaut=max(max(abs(wautnew-waut)));
    waut=wautnew;
end
for j=1:state
    for g=1:agent
    if waut(j,g)>0
waut(j,g)=waut(j,g)*(1-penalty);
    else
waut(j,g)=waut(j,g)/(1-penalty);
    end
    end
end

% Solve problem without enforcement constraints at every grid point

w=zeros(n,state,agent);
gammar=zeros(n,state,agent);
c=zeros(n,state,agent);

% Note that phi contains the Lagrange multipliers on the enforcement
% constraint
phi=zeros(n,state,agent);

count=0

for j=1:state
    for g=1:agent
c(:,j,g)=Y(j)*gammagrid(:,agent+g)./(1+gammagrid(:,4));
     end
end
cfirstbest=c;
csolold=cfirstbest;
gammar=zeros(n,state,agent);
for j=1:state
for k=1:agent   
gammar(:,j,k)=gammagrid(:,agent+k);

end
end
% Compute the value functions
w=zeros(n,state,agent);
wnew=zeros(n,state,agent);
dw=1;
while dw>crit
    for j=1:state
        wnew(:,j,1)=utility(rho,cfirstbest(:,j,1))+beta*w(:,:,1)*P(j,:)';
        wnew(:,j,2)=utility(rho,cfirstbest(:,j,2))+beta*w(:,:,2)*P(j,:)';
        

    end
            dw=max(max(max(abs(w-wnew))));
        w=wnew;
end

wfirstbest=w;
%%%
ewfirstbest=zeros(n,state,agent);
for j=1:state
        for g=1:agent
ewfirstbest(:,j,g)=wfirstbest(:,:,g)*P(j,:)';
        end
end
ewsb=ewfirstbest;
wsb=wfirstbest;
wsbnew=wsb;


gap=1;
gapend=0.001;
options=optimset('Display','iter');
scalefactor=1;

count=0;
eps=0.0000001;
t=0;
csol=[]; wsbnewsol=[]; phisol=[];
while gap>0.001
for j=1:state
   
%%%%Constraint binding for person 1
rr=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)<waut(j,1) & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2)); 
if size(rr)>0
index=[rr(end) j];       
      
[xsol,output]=agent1solverho_2(index,ewsb,waut);

csol(rr,j,1)=xsol(1);
csol(rr,j,2)=Y(j)-xsol(1);
gammar(rr,j,2)=csol(rr,j,2)./csol(rr,j,1);
for g=1:agent
wsbnewsol(rr,j,g)=output(g);
end
phisol(rr,j,1)=(gammagrid(rr,4)./gammar(rr,j,2)).^rho-1;
phisol(rr,j,2)=0;
end
%%%%Constraint binding for person 2
rr=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)>=waut(j,1) & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)<waut(j,2)); 
if size(rr)>0
index=[rr(1) j];
[xsol,output]=agent2solverho_2(index,ewsb,waut);

csol(rr,j,1)=xsol(1);
csol(rr,j,2)=Y(j)-xsol(1);
gammar(rr,j,2)=csol(rr,j,2)./csol(rr,j,1);
for g=1:agent
wsbnewsol(rr,j,g)=output(g);
end
phisol(rr,j,2)=1/(gammagrid(rr,4)./gammar(rr,j,2)).^rho-1;
phisol(rr,j,2)=(gammar(rr,j,2)./gammagrid(rr,4)).^rho-1;
phisol(rr,j,1)=0;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Constraint binding for noone%%%%%%%
rr=find( utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)>=waut(j,1) & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2) )
if size(rr)>0        

wsbnewsol(rr,j,1)=utility(rho,cfirstbest(rr,j,1))+beta*ewsb(rr,j,1);
wsbnewsol(rr,j,2)=utility(rho,cfirstbest(rr,j,2))+beta*ewsb(rr,j,2);

csol(rr,j,1)=cfirstbest(rr,j,1);
csol(rr,j,2)=cfirstbest(rr,j,2);
gammar(rr,j,2)=csol(rr,j,2)./csol(rr,j,1);
phisol(rr,j,1)=0;
phisol(rr,j,2)=0;

end

end % end of j loop


for g=1:agent
gap(g)=scalefactor*norm((wsbnewsol(:,:,g)-wsb(:,:,g)));
end
for j=1:state
for g=1:agent
gammar(:,j,g)=csol(:,j,g)./csol(:,j,1);
end
end
gap=sum(gap)
gap1=max(max(max(abs(invutility(beta,rho,wsbnewsol(:,:,1))-invutility(beta,rho,wsb(:,:,1)))./invutility(beta,rho,wsb(:,:,1)))),max(max(abs(invutility(beta,rho,wsbnewsol(:,:,2))-invutility(beta,rho,wsb(:,:,2)))./invutility(beta,rho,wsb(:,:,2)))))/(1-beta);
gap2=max(max(max(abs(csol(:,:,1)-csolold(:,:,1))./csol(:,:,1))),max(max(abs(csol(:,:,2)-csolold(:,:,2))./csol(:,:,2))));
gap=max(gap1,gap2)
wsb=wsbnewsol;
   
for j=1:state
    for g=1:agent
ewsb(:,j,g)=wsb(:,:,g)*P(j,:)';
    end
end
csolold=csol;

end
PROBDEV2(1:state,1:state,1)=P;
STATEDEV2(1:state,1:agent,1)=S;
WDEV2(1:size(wsb,1),1:size(wsb,2),1:2,Qbeta,Qrho)=wsb;
EWDEV2(1:size(ewsb,1),1:size(ewsb,2),1:2,Qbeta,Qrho)=ewsb;
PHIDEV2(1:size(phisol,1),1:size(phisol,2),1:2,Qbeta,Qrho)=phisol;
CDEV2(1:size(csol,1),1:size(csol,2),1:2,Qbeta,Qrho)=csol;



save ('RESULTS2ALL.mat', 'PROBDEV2','STATEDEV2','WDEV2','EWDEV2','PHIDEV2','CDEV2');


